/* Copyright 2007-2009 dnAnalytics Project.
 *
 * Contributors to this file:
 * Marcus Cuda
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
#include <algorithm>
#include "mkl_cblas.h"
#include "mkl_lapack.h"
#include "mkl_vml.h"
#include "mkl_trans.h"
#include "linear_algebra.h"
using namespace std;

extern "C" {

	size_t s_index_max(const int n, const float data[]){
		return cblas_isamax( n, data, 1);
	}

	size_t d_index_max(const int n, const double data[]){
		return cblas_idamax( n, data, 1);
	}

	size_t c_index_max(const int n, const MKL_Complex8 data[]){
		return cblas_icamax( n, data, 1);
	}

	size_t z_index_max(const int n, const MKL_Complex16 data[]){
		return cblas_izamax( n, data, 1);
	}

	size_t s_index_min(const int n, const float data[]){
		return cblas_isamin( n, data, 1);
	}

	size_t d_index_min(const int n, const double data[]){
		return cblas_idamin( n, data, 1);
	}

	size_t c_index_min(const int n, const MKL_Complex8 data[]){
		return cblas_icamin( n, data, 1);
	}

	size_t z_index_min(const int n, const MKL_Complex16 data[]){
		return cblas_izamin( n, data, 1);
	}

	float s_absolute_sum(const int n, const float data[]){
		return cblas_sasum(n, data, 1);
	}

	double d_absolute_sum(const int n, const double data[]){
		return cblas_dasum(n, data, 1);
	}

	float c_absolute_sum(const int n, const MKL_Complex8 data[]){
		return cblas_scasum(n, data, 1);
	}

	double z_absolute_sum(const int n, const MKL_Complex16 data[]){
		return cblas_dzasum(n, data, 1);
	}

	float s_dot_product(const int n, const float x[], const float y[]){
		return cblas_sdot(n, x, 1, y, 1);
	}

	double d_dot_product(const int n, const double x[], const double y[]){
		return cblas_ddot(n, x, 1, y, 1);
	}

	MKL_Complex8 c_dot_product(const int n, const MKL_Complex8 x[], const MKL_Complex8 y[]){
		MKL_Complex8 ret;
		cblas_cdotu_sub(n, x, 1, y, 1, &ret);
		return ret;
	}

	MKL_Complex16 z_dot_product(const int n, const MKL_Complex16 x[], const MKL_Complex16 y[]){
		MKL_Complex16 ret;
		cblas_zdotu_sub(n, x, 1, y, 1, &ret);
		return ret;
	}

	MKL_Complex8 c_dot_product_c(const int n, const MKL_Complex8 x[], const MKL_Complex8 y[]){
		MKL_Complex8 ret;
		cblas_cdotc_sub(n, x, 1, y, 1, &ret);
		return ret;
	}

	MKL_Complex16 z_dot_product_c(const int n, const MKL_Complex16 x[], const MKL_Complex16 y[]){
		MKL_Complex16 ret;
		cblas_zdotu_sub(n, x, 1, y, 1, &ret);
		return ret;
	}

	void s_scale(const int n, const float alpha, float x[]){
		cblas_sscal(n, alpha, x, 1);
	}

	void d_scale(const int n, const double alpha, double x[]){
		cblas_dscal(n, alpha, x, 1);
	}

	void c_scale(const int n, const MKL_Complex8 alpha, MKL_Complex8 x[]){
		cblas_cscal(n, &alpha, x, 1);
	}
	
	void z_scale(const int n, const MKL_Complex16 alpha, MKL_Complex16 x[]){
		cblas_zscal(n, &alpha, x, 1);
	}
	
	void s_axpy( const int n, const float alpha, const float x[], float y[]){
		cblas_saxpy(n, alpha, x, 1, y, 1);
	}

	void d_axpy( const int n, const double alpha, const double x[], double y[]){
		cblas_daxpy(n, alpha, x, 1, y, 1);
	}

	void c_axpy( const int n, const MKL_Complex8 alpha, const MKL_Complex8 x[], MKL_Complex8 y[]){
		cblas_caxpy(n, &alpha, x, 1, y, 1);
	}

	void z_axpy( const int n, const MKL_Complex16 alpha, const MKL_Complex16 x[], MKL_Complex16 y[]){
		cblas_zaxpy(n, &alpha, x, 1, y, 1);
	}

	void d_vector_add( const int n, const double x[], const double y[], double ret[]){
		vdAdd( n, x, y, ret );
	}

	void d_vector_subtract( const int n, const double x[], const double y[], double ret[]){
		vdSub( n, x, y, ret );
	}

	void d_vector_multiply( const int n, const double x[], const double y[], double ret[]){
		vdMul( n, x, y, ret );
	}

	void s_matrix_multiply(const int m, const int n, const int k, const float x[],
		const float y[], float c[]){
			cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1.0, x, m, y, k, 0.0, c, m);
	}

	void d_matrix_multiply(const int m, const int n, const int k, const double x[],
		const double y[], double c[]){
			cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1.0, x, m, y, k, 0.0, c, m);
	}

	void c_matrix_multiply(const int m, const int n, const int k, const MKL_Complex8 x[],
		const MKL_Complex8 y[], MKL_Complex8 c[]){
		MKL_Complex8 one;
	    one.real = 1;
		MKL_Complex8 zero;
	
		cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, &one, x, m, y, k, &zero, c, m);
	}

	void z_matrix_multiply(const int m, const int n, const int k, const MKL_Complex16 x[],
		const MKL_Complex16 y[], MKL_Complex16 c[]){
		MKL_Complex16 one;
	    one.real = 1;
		MKL_Complex16 zero;	
		cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, &one, x, m, y, k, &zero, c, m);
	}
	
	int s_lu_factor(int m, float a[], int ipiv[]){
		int info;
		SGETRF(&m,&m,a,&m,ipiv,&info);
		for(int i = 0; i < m; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	int d_lu_factor(int m, double a[], int ipiv[]){
		int info;

		DGETRF(&m,&m,a,&m,ipiv,&info);
		for(int i = 0; i < m; ++i ){
			ipiv[i] -= 1;
		}

		return info;
	}

	int c_lu_factor(int m, MKL_Complex8 a[], int ipiv[]){
		int info;
		CGETRF(&m,&m,a,&m,ipiv,&info);
		for(int i = 0; i < m; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	int z_lu_factor(int m, MKL_Complex16 a[], int ipiv[]){
		int info;
		ZGETRF(&m,&m,a,&m,ipiv,&info);
		for(int i = 0; i < m; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	int s_lu_inverse( int n, float a[], int ipiv[] ){
		int i;
		for(i= 0; i < n; ++i ){
			ipiv[i] += 1;
		}
		int info;
		int ispec = 1;
		int n2 = -1;
		int bs = ILAENV(&ispec, "SGETRI", "", &n, &n2, &n2, &n2); 
		int lwork = bs*n;
		float* work = new float[lwork];

		SGETRI(&n,a,&n,ipiv,work,&lwork,&info);

		delete[] work;
		for(i = 0; i < n; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	int d_lu_inverse( int n, double a[], int ipiv[] ){
		int i;
		for(i= 0; i < n; ++i ){
			ipiv[i] += 1;
		}
		int info;
		int ispec = 1;
		int n2 = -1;
		int bs = ILAENV(&ispec, "DGETRI", "", &n, &n2, &n2, &n2); 
		int lwork = bs*n;
		double* work = new double[lwork];

		DGETRI(&n,a,&n,ipiv,work,&lwork,&info);

		delete[] work;
		for(i = 0; i < n; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	int c_lu_inverse( int n, MKL_Complex8 a[], int ipiv[] ){
		int i;
		for(i= 0; i < n; ++i ){
			ipiv[i] += 1;
		}
		int info;
		int ispec = 1;
		int n2 = -1;
		int bs = ILAENV(&ispec, "CGETRI", "", &n, &n2, &n2, &n2); 
		int lwork = bs*n;
		MKL_Complex8* work = new MKL_Complex8[lwork];

		CGETRI(&n,a,&n,ipiv,work,&lwork,&info);

		delete[] work;
		for(i = 0; i < n; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	int z_lu_inverse( int n, MKL_Complex16 a[], int ipiv[] ){
		int i;
		for(i= 0; i < n; ++i ){
			ipiv[i] += 1;
		}
		int info;
		int ispec = 1;
		int n2 = -1;
		int bs = ILAENV(&ispec, "ZGETRI", "", &n, &n2, &n2, &n2); 
		int lwork = bs*n;
		MKL_Complex16* work = new MKL_Complex16[lwork];

		ZGETRI(&n,a,&n,ipiv,work,&lwork,&info);

		delete[] work;
		for(i = 0; i < n; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	int s_lu_solve(int n, int nrhs, float a[], int ipiv[], float b[]){
		int info;
		int i;    
		for(i = 0; i < n; ++i ){
			ipiv[i] += 1;
		}

		char trans = 'N';
		SGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info);
		for(i = 0; i < n; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	int d_lu_solve(int n, int nrhs, double a[], int ipiv[], double b[]){
		int info;
		int i;    
		for(i = 0; i < n; ++i ){
			ipiv[i] += 1;
		}

		char trans = 'N';
		DGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info);
		for(i = 0; i < n; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	int c_lu_solve(int n, int nrhs, MKL_Complex8 a[], int ipiv[], MKL_Complex8 b[]){
		int info;
		int i;    
		for(i = 0; i < n; ++i ){
			ipiv[i] += 1;
		}

		char trans = 'N';
		CGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info);
		for(i = 0; i < n; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	int z_lu_solve(int n, int nrhs, MKL_Complex16 a[], int ipiv[], MKL_Complex16 b[]){
		int info;
		int i;    
		for(i = 0; i < n; ++i ){
			ipiv[i] += 1;
		}

		char trans = 'N';
		ZGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info);
		for(i = 0; i < n; ++i ){
			ipiv[i] -= 1;
		}
		return info;
	}

	float s_matrix_norm(char norm, int m, int n, float a[])
	{
		float* work = new float[m];
		float ret = SLANGE(&norm, &m, &n, a, &m, work);
		delete[] work;
		return ret;
	}

	double d_matrix_norm(char norm, int m, int n, double a[])
	{
		double* work = new double[m];
		double ret = DLANGE(&norm, &m, &n, a, &m, work);
		delete[] work;
		return ret;
	}
	
	float c_matrix_norm(char norm, int m, int n, MKL_Complex8 a[])
	{
		float* work = new float[m];
		float ret = CLANGE(&norm, &m, &n, a, &m, work);
		delete[] work;
		return ret;
	}
	
	double z_matrix_norm(char norm, int m, int n, MKL_Complex16 a[])
	{
		double* work = new double[m];
		double ret = ZLANGE(&norm, &m, &n, a, &m, work);
		delete[] work;
		return ret;
	}

	int s_condition(char type, int n, float norm, float a[], float* rcond){
		int info = 0;
		float* work = new float[4 * n];
		int* iwork = new int[n];
		SGECON(&type,&n,a,&n,&norm, rcond, work, iwork, &info);
		delete[] work;
		delete[] iwork;
		return info;
	}
	
	int d_condition(char type, int n, double norm, double a[], double* rcond){
		int info = 0;
		double* work = new double[4 * n];
		int* iwork = new int[n];
		DGECON(&type,&n,a,&n,&norm, rcond, work, iwork, &info);
		delete[] work;
		delete[] iwork;
		return info;
	}

	int c_condition(char type, int n, float norm, MKL_Complex8 a[], float* rcond){
		int info = 0;
		MKL_Complex8* work = new MKL_Complex8[4 * n];
		float* rwork = new float[2*n];
		CGECON(&type,&n,a,&n,&norm, rcond, work, rwork, &info);
		delete[] work;
		delete[] rwork;
		return info;
	}

	int z_condition(char type, int n, double norm, MKL_Complex16 a[], double* rcond){
		int info = 0;
		MKL_Complex16* work = new MKL_Complex16[4 * n];
		double* rwork = new double[2*n];
		ZGECON(&type,&n,a,&n,&norm, rcond, work, rwork, &info);
		delete[] work;
		delete[] rwork;
		return info;
	}

	int s_cholesky_factor(int n, float a[]){
		char uplo = 'L';
		int info = 0;
		SPOTRF(&uplo, &n, a, &n, &info);
		for (int i = 0; i < n; ++i)
		{
			int index = i * n;
			for (int j = 0; j < n && i > j; ++j)
			{
				a[index + j] = 0;
			}
		}
		return info;
	}

	int d_cholesky_factor(int n, double* a){
		char uplo = 'L';
		int info = 0;
		DPOTRF(&uplo, &n, a, &n, &info);
		for (int i = 0; i < n; ++i)
		{
			int index = i * n;
			for (int j = 0; j < n && i > j; ++j)
			{
				a[index + j] = 0;
			}
		}
		return info;
	}

	int c_cholesky_factor(int n, MKL_Complex8 a[]){
		char uplo = 'L';
		int info = 0;
		MKL_Complex8 zero;
		zero.real = 0.0;
		zero.real = 0.0;
		CPOTRF(&uplo, &n, a, &n, &info);
		for (int i = 0; i < n; ++i)
		{
			int index = i * n;
			for (int j = 0; j < n && i > j; ++j)
			{
				a[index + j] = zero;
			}
		}
		return info;
	}

	int z_cholesky_factor(int n, MKL_Complex16 a[]){
		char uplo = 'L';
		int info = 0;
		MKL_Complex16 zero;
		zero.real = 0.0;
		zero.real = 0.0;
		ZPOTRF(&uplo, &n, a, &n, &info);
		for (int i = 0; i < n; ++i)
		{
			int index = i * n;
			for (int j = 0; j < n && i > j; ++j)
			{
				a[index + j] = zero;
			}
		}
		return info;
	}

	int s_cholesky_solve(int n, int nrhs, float a[], float b[]){
		char uplo = 'L';
		int info = 0;
		SPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info);
		return info;
	}

	int d_cholesky_solve(int n, int nrhs, double a[], double b[]){
		char uplo = 'L';
		int info = 0;
		DPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info);
		return info;
	}

	int c_cholesky_solve(int n, int nrhs, MKL_Complex8 a[], MKL_Complex8 b[]){
		char uplo = 'L';
		int info = 0;
		CPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info);
		return info;
	}

	int z_cholesky_solve(int n, int nrhs, MKL_Complex16 a[], MKL_Complex16 b[]){
		char uplo = 'L';
		int info = 0;
		ZPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info);
		return info;
	}

	int s_qr_factor(int m, int n, float r[], float tau[], float q[]){
		int minus_one = -1;
		int ispec = 1;
		int bs = ILAENV(&ispec, "SGEQRF", "", &m, &n, &minus_one, &minus_one);
		int len = n * bs;
		float* work = new float[len];
		int info = 0;
		SGEQRF(&m, &n, r, &m, tau, work, &len, &info);
		delete[] work;

		for (int i = 0; i < m; ++i)
		{
			for (int j = 0; j < m && j < n; ++j)
			{
				if (i > j)
				{
					q[j * m + i] = r[j * m + i];
				}
			}
		}

		//compute the q elements explicitly
		if (m <= n)
		{
			bs = ILAENV(&ispec, "SORGQR", "", &m, &m, &m, &minus_one);
			len = n * bs;
			work = new float[len];
			SORGQR(&m, &m, &m, q, &m, tau, work, &len, &info);
			delete[] work;
		}
		else
		{
			bs = ILAENV(&ispec, "SORGQR", "", &m, &m, &n, &minus_one);
			len = n * bs;
			work = new float[len];
			SORGQR(&m, &n, &n, q, &m, tau, work, &len, &info);
			delete[] work;
		}
		return info;
	}

	int d_qr_factor(int m, int n, double r[], double tau[], double q[]){
		int minus_one = -1;
		int ispec = 1;
		int bs = ILAENV(&ispec, "DGEQRF", "", &m, &n, &minus_one, &minus_one);
		int len = n * bs;
		double* work = new double[len];
		int info = 0;
		DGEQRF(&m, &n, r, &m, tau, work, &len, &info);
		delete[] work;

		for (int i = 0; i < m; ++i)
		{
			for (int j = 0; j < m && j < n; ++j)
			{
				if (i > j)
				{
					q[j * m + i] = r[j * m + i];
				}
			}
		}

		//compute the q elements explicitly
		if (m <= n)
		{
			bs = ILAENV(&ispec, "DORGQR", "", &m, &m, &m, &minus_one);
			len = n * bs;
			work = new double[len];
			DORGQR(&m, &m, &m, q, &m, tau, work, &len, &info);

			delete[] work;
		}
		else
		{
			bs = ILAENV(&ispec, "DORGQR", "", &m, &m, &n, &minus_one);

			len = n * bs;
			work = new double[len];
			DORGQR(&m, &n, &n, q, &m, tau, work, &len, &info);
			delete[] work;
		}
		return info;
	}

	int c_qr_factor(int m, int n, MKL_Complex8 r[], MKL_Complex8 tau[], MKL_Complex8 q[]){
		int minus_one = -1;
		int ispec = 1;
		int bs = ILAENV(&ispec, "CGEQRF", "", &m, &n, &minus_one, &minus_one);
		int len = n * bs;
		MKL_Complex8* work = new MKL_Complex8[len];
		int info = 0;
		CGEQRF(&m, &n, r, &m, tau, work, &len, &info);
		delete[] work;

		for (int i = 0; i < m; ++i)
		{
			for (int j = 0; j < m && j < n; ++j)
			{
				if (i > j)
				{
					q[j * m + i] = r[j * m + i];
				}
			}
		}

		//compute the q elements explicitly
		if (m <= n)
		{
			bs = ILAENV(&ispec, "CUNGQR", "", &m, &m, &m, &minus_one);
			len = n * bs;
			work = new MKL_Complex8[len];
			CUNGQR(&m, &m, &m, q, &m, tau, work, &len, &info);
			delete[] work;
		}
		else
		{
			bs = ILAENV(&ispec, "CUNGQR", "", &m, &m, &n, &minus_one);
			len = n * bs;
			work = new MKL_Complex8[len];
			CUNGQR(&m, &n, &n, q, &m, tau, work, &len, &info);
			delete[] work;
		}
		return info;
	}

	int z_qr_factor(int m, int n, MKL_Complex16 r[], MKL_Complex16 tau[], MKL_Complex16 q[]){
		int minus_one = -1;
		int ispec = 1;
		int bs = ILAENV(&ispec, "ZGEQRF", "", &m, &n, &minus_one, &minus_one);
		int len = n * bs;
		MKL_Complex16* work = new MKL_Complex16[len];
		int info = 0;
		ZGEQRF(&m, &n, r, &m, tau, work, &len, &info);
		delete[] work;

		for (int i = 0; i < m; ++i)
		{
			for (int j = 0; j < m && j < n; ++j)
			{
				if (i > j)
				{
					q[j * m + i] = r[j * m + i];
				}
			}
		}

		//compute the q elements explicitly
		if (m <= n)
		{
			bs = ILAENV(&ispec, "ZUNGQR", "", &m, &m, &m, &minus_one);
			len = n * bs;
			work = new MKL_Complex16[len];
			ZUNGQR(&m, &m, &m, q, &m, tau, work, &len, &info);
			delete[] work;
		}
		else
		{
			bs = ILAENV(&ispec, "ZUNGQR", "", &m, &m, &n, &minus_one);
			len = n * bs;
			work = new MKL_Complex16[len];
			ZUNGQR(&m, &n, &n, q, &m, tau, work, &len, &info);
			delete[] work;
		}
		return info;
	}

	int s_qr_solve(int m, int n, int bn, float r[], float b[], float tau[], float x[]){
		int minus_one = -1;
		int ispec = 1;
		int bs = ILAENV(&ispec, "SORMQR", "LT", &m, &bn, &n, &minus_one);
		int len = n * bs;
		float* work = new float[len];
		char side ='L';
		char tran = 'T';
		int info = 0;
		SORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info);
		delete[] work;
		cblas_strsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, 1.0, r, m, b, m);
		for (int i = 0; i < n; ++i)
		{
			for (int j = 0; j < bn; ++j)
			{
				x[j * n + i] = b[j * m + i];
			}
		}
		return info;
	}

	int d_qr_solve(int m, int n, int bn, double r[], double b[], double tau[], double x[]){
		int minus_one = -1;
		int ispec = 1;
		int bs = ILAENV(&ispec, "DORMQR", "LT", &m, &bn, &n, &minus_one);
		int len = n * bs;
		double* work = new double[len];
		char side ='L';
		char tran = 'T';
		int info = 0;
		DORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info);
		delete[] work;
		cblas_dtrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, 1.0, r, m, b, m);
		for (int i = 0; i < n; ++i)
		{
			for (int j = 0; j < bn; ++j)
			{
				x[j * n + i] = b[j * m + i];
			}
		}
		return info;
	}

	int c_qr_solve(int m, int n, int bn, MKL_Complex8 r[], MKL_Complex8 b[], MKL_Complex8 tau[], MKL_Complex8 x[]){
		int minus_one = -1;
		int ispec = 1;
		int bs = ILAENV(&ispec, "CUNMQR", "LT", &m, &bn, &n, &minus_one);
		int len = n * bs;
		MKL_Complex8* work = new MKL_Complex8[len];
		char side ='L';
		char tran = 'T';
		int info = 0;
		CUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info);
		delete[] work;
		MKL_Complex8 one;
		one.real = 1.0;
		cblas_ctrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, &one, r, m, b, m);
		for (int i = 0; i < n; ++i)
		{
			for (int j = 0; j < bn; ++j)
			{
				x[j * n + i] = b[j * m + i];
			}
		}
		return info;
	}

	int z_qr_solve(int m, int n, int bn, MKL_Complex16 r[], MKL_Complex16 b[], MKL_Complex16 tau[], MKL_Complex16 x[]){
		int minus_one = -1;
		int ispec = 1;
		int bs = ILAENV(&ispec, "ZUNMQR", "LT", &m, &bn, &n, &minus_one);
		int len = n * bs;
		MKL_Complex16* work = new MKL_Complex16[len];
		char side ='L';
		char tran = 'T';
		int info = 0;
		ZUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info);
		delete[] work;
		MKL_Complex16 one;
		one.real = 1.0;
		cblas_ztrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, &one, r, m, b, m);
		for (int i = 0; i < n; ++i)
		{
			for (int j = 0; j < bn; ++j)
			{
				x[j * n + i] = b[j * m + i];
			}
		}
		return info;
	}

	int s_svd_factor(bool compute_vectors, int m, int n, float a[], float s[], float u[], float v[]){
		int len = max((3 * min(m, n) + max(m, n)), 5 * min(m, n));
		float* work = new float[len];
		int info = 0;
		char job = compute_vectors ? 'A' : 'N';
		SGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info);
		delete[] work;
		return info;
	}
	
	int d_svd_factor(bool compute_vectors, int m, int n, double a[], double s[], double u[], double v[]){
		int len = max((3 * min(m, n) + max(m, n)), 5 * min(m, n));
		double* work = new double[len];
		int info = 0;
		char job = compute_vectors ? 'A' : 'N';
		DGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info);
		delete[] work;
		return info;
	}

	int c_svd_factor(bool compute_vectors, int m, int n, MKL_Complex8 a[], float s[], MKL_Complex8 u[], MKL_Complex8 v[]){
		int len = max((3 * min(m, n) + max(m, n)), 5 * min(m, n));
		MKL_Complex8* work = new MKL_Complex8[len];
		int info = 0;
		char job = compute_vectors ? 'A' : 'N';
		float* rwork = new float[5*min(m,n)];
		CGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, rwork, &info);
		delete[] work;
		delete[] rwork;
		return info;
	}
	
	int z_svd_factor(bool compute_vectors, int m, int n, MKL_Complex16 a[], double s[], MKL_Complex16 u[], MKL_Complex16 v[]){
		int len = max((3 * min(m, n) + max(m, n)), 5 * min(m, n));
		MKL_Complex16* work = new MKL_Complex16[len];
		int info = 0;
		char job = compute_vectors ? 'A' : 'N';
		double* rwork = new double[5*min(m,n)];
		ZGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, rwork, &info);
		delete[] work;
		return info;
	}

	double d_householder(int n, double x[]){
		double tau = 0;
		int incx = 1;
		DLARFG(&n, x, x+1, &incx, &tau);

		return tau;
	}

	void d_transpose(int rows, int cols, double a[], double b[]){
		mkl_domatcopy('C', 'T', rows, cols, 1.0, a, rows, b, cols);
	}

	void d_gemm(const int m, const int n, const int k, int transA, int transB, const double alpha, const double x[], const double y[], const double beta, double c[]){
		
		CBLAS_TRANSPOSE transposeA = transA ? CblasTrans : CblasNoTrans;
		CBLAS_TRANSPOSE transposeB = transB ? CblasTrans : CblasNoTrans;
		int lda = transA ? k : m;
		int ldb = transB ? n : k;
		
		cblas_dgemm(CblasColMajor, transposeA, transposeB, m, n, k, alpha, x, lda, y, ldb, beta, c, m);
	}
}
